'''
Glue for returning descriptive statistics.
'''
import numpy as np
from scipy import stats

# need the dir containing scipystats in your PYTHONPATH
# file structure WILL change
from scipystats.sandbox.load_dataset import load_dataset
from scipystats.sandbox.string2dummy import string2dummy as s2d


#############################################
#
#============================================
#       Univariate Descriptive Statistics
#============================================
#

def sign_test(samp,Mu0=0):
        '''
        This is just a simple signs test with Mu0=0 by default (though
        the median is often used in practice)
        Not sure if it is already in numpy or scipy.
        I couldn't find it.

        Returns M, p-value
        where
        M=(N(+) - N(-))/2, N(+) is the number of values above Mu0,
        N(-) is the number of values below.  Values equal to Mu0 
        are discarded.

        The p-value for M is calculated using the binomial distrubution
        and can be intrepreted the same as for a t-test.

        See Also
        ---------
        scipy.stats.wilcoxon

        Notes
        ------
        Needs to be reviewed and tested.
        '''
        pos=np.sum(samp>Mu0)
        neg=np.sum(samp<Mu0)
        M=(pos-neg)/2.
        p=stats.binom_test(min(pos,neg),pos+neg,.5)
        return M,p


def descstats(x,v):
    '''
    Prints descriptive statistics for one or multiple variables.

    Parameters
    ------------
    x: numpy array
        `x` is the array that holds all of your data.
        It can be any type of array.
    v: list
        `v` is a list of the field names of variables.        

    Example
    ----------
    decstats(data_array,['x_1','x_2','x_3'])

    '''
#   todo: check that x and v are not empty
#   todo: check that var in v exists
#   todo: scipy.stats.mode? (http://projects.scipy.org/scipy/ticket/905)
    # don't worry about premature optimization...
#   todo: add optional user-defined precision?
#   todo: throw different numbers at it to test string formatting
#   todo: should this be a class that prints, but also holds data?
#   todo: have a closer look at matplotlib.mlab.group_by for stats by attribute
#   todo: explicit handling of missing values (this has to be handled consistently
#         at the data loading level)
#   use scipy.stats.describe?
        
    if (len(v)==1):
        print '''
    ---------------------------------------------
    Univariate Descriptive Statistics
    ---------------------------------------------

    Var. Name   %(name)12s
    ----------
    Obs.          %(nobs)22i  Range              %(range)22s
    Sum of Wts.   %(sum)22s  Coeff. of Variation     %(coeffvar)22.4g
    Mode          %(mode)22.4g  Skewness                %(skewness)22.4g
    Repeats       %(nmode)22i  Kurtosis                %(kurtosis)22.4g
    Mean          %(mean)22.4g  Uncorrected SS          %(uss)22.4g
    Median        %(median)22.4g  Corrected SS            %(ss)22.4g
    Variance      %(variance)22.4g  Sum Observations        %(sobs)22.4g   
    Std. Dev.     %(stddev)22.4g
    ''' % {'name': v[0], 'sum': 'N/A', 'nobs': len(x[v[0]]), 'mode': \
    stats.mode(x[v[0]])[0][0], 'nmode': stats.mode(x[v[0]])[1][0], \
    'mean': x[v[0]].mean(), 'median': np.median(x[v[0]]), 'range': \
    '('+str(x[v[0]].min())+', '+str(x[v[0]].max())+')', 'variance': \
    x[v[0]].var(), 'stddev': x[v[0]].std(), 'coeffvar': \
    stats.variation(x[v[0]]), 'skewness': stats.skew(x[v[0]]), \
    'kurtosis': stats.kurtosis(x[v[0]]), 'uss': stats.ss(x[v[0]]),\
    'ss': stats.ss(x[v[0]]-x[v[0]].mean()), 'sobs': np.sum(x[v[0]])}

        print '''
    Percentiles
    -------------
    1  %%          %12.4g
    5  %%          %12.4g 
    10 %%          %12.4g
    25 %%          %12.4g

    50 %%          %12.4g

    75 %%          %12.4g
    90 %%          %12.4g
    95 %%          %12.4g
    99 %%          %12.4g
    ''' % tuple([stats.scoreatpercentile(x[v[0]],per) for per in (1,5,10,25,50,75,90,95,99)])
        t,p_t=stats.ttest_1samp(x[v[0]],0) 
        M,p_M=sign_test(x[v[0]])
        S,p_S=stats.wilcoxon(x[v[0]])

        print '''

    Tests of Location (H0: Mu0=0)
    -----------------------------
    Test                Statistic       Two-tailed probability
    -----------------+-----------------------------------------
    Student's t      |  t %7.5f   Pr > |t|   <%.4f
    Sign             |  M %8.1f   Pr >= |M|  <%.4f
    Signed Rank      |  S %8.1f   Pr >= |S|  <%.4f

    ''' % (t,p_t,M,p_M,S,p_S)
# Should this be part of a 'descstats'
# in any event these should be split up, so that they can be called
# individually and only returned together if someone calls summary
# or something of the sort
    
    elif (len(v)>1):  
        summary='''
    Var. Name   |     Obs.        Mean    Std. Dev.           Range
    ------------+--------------------------------------------------------'''
        for vr in v: 
            print "%(name)15s %(obs)9i %(mean)12.4g %(stddev)12.4g %(range)20s" %\
            {'name': vr, 'obs': len(x[vr]), 'mean': x[vr].mean(), 'stddev': \
            x[vr].std(), 'range': '('+str(x[vr].min())+', '+str(x[vr].max())+')'} 

if __name__=='__main__':
# test descstats
    import os
    loc='http://eagle1.american.edu/~js2796a/data/handguns_data.csv'
    relpath=(load_dataset(loc))
    dta=np.recfromcsv(relpath)  
#    descstats(dta,['stpop'])
#    raw_input('Hit enter for multivariate test')
#    descstats(dta,['stpop','avginc','vio'])

# with plain arrays
    import string2dummy as s2d
    dts=s2d.string2dummy(dta)
    ndts=np.vstack(dts[col] for col in dts.dtype.names)
# observations in columns and data in rows
# is easier for the call to stats
    
# what to make of  
# ndts=np.column_stack(dts[col] for col in dts.dtype.names)
# ntda=ntds.swapaxis(1,0)
# ntda is ntds returns false?


# What about the stats for the string variables, that whole array should 
# be a "dummy array"

# or now we just have detailed information about the different strings
# would this approach ever be inappropriate for a string typed variable
# other than dates?
    descstats(ndts, [1])
    raw_input("Enter to try second part")
    descstats(ndts, [1,20,3])
